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Abstract 

A  diffiisioo  approximation  for  a  network  of  continuous  time  reservoirs  with  power  law  release  rules  is  examined. 
Under  a  mild  assumption  on  the  inflow  processes,  we  show  that  for  physically  reasonable  \’alues  of  the  power  law 
constants,  the  system  of  processes  converges  to  a  multi-dimensional  Gaussiaa  diffusion  process.  We  also  illustrate  how 
the  limiting  Gaussian  process  may  be  used  to  compute  approximations  to  the  original  system  of  processes.  In  addition, 
we  study  the  quality  of  our  approximations  by  comparing  them  to  results  obtained  by  simulations  of  the  original  watershed 
model.  The  simulations  offer  support  for  the  use  of  the  approximatioo  developed  here. 


1.  Introduction 

The  linear  reservoir  plays  a  central  role  in  hydrology  (Nash,  1957, 1959).  Although  this 
model  has  been  successfully  used  for  many  applications,  a  more  general  system  of  reservoirs 
with  a  power  law  storage-runoff  relation  can  more  accurately  represent  the  physical 
characteristics  of  a  natural  storage  system.  The  nonlinear  reservoir  has  been  examined  by 
Laurensen  (1964),  Mein  et  al  (1974),  K]eme$  (1978),  Hughes  and  Murrel  (1986),  and  Glynn 
(1989))  among  others.  The  natural  extension  of  this  model  to  a  nonlinear  cascade  of  reservoirs 
has  also  been  investigated  (KlemeS  et  al.  (1975,1985)).  Unfortunately  closed  form  formulae  for 
the  statistical  properties  of  these  reservoirs  are  only  available  for  special  cases  (e.g.  the  linear 
reservoir,  the  fixed  release  reservoir,  the  power  law  release  with  special  inflow  distribution  and 
exponent  less  than  one). 

Two  approaches  have  generally  been  used  to  model  the  continuous  time  reservoir.  In  the 
first  case,  the  inflows  are  directly  modeled  as  diffusion  processes  and  the  problem  requires 
solving  nonlinear  stochastic  differential  equations  (Unny  (1984),  Unny  and  Karmeshu  (1984)). 
Such  models  have  also  been  used  to  study  optimal  control  problems  associated  with  reservoirs 
(Harrison  and  Taylor  (1978),  Pliska  (1975)).  The  second  approach  assumes  that  inflows  follow 
a  compound  Poisson  process.  The  linear  cascade  with  general  additive  homogeneous  inflows  has 
been  investigated  by  Moran  (1967).  Exact  solutions  for  the  single  reservoir  with  a  power  law 
release  rule  are  available  for  the  case  of  compound  Poisson  inflows  where  the  power  is  less  than 
one  (Harrison  and  Resnidc  (1976)).  Approximate  solutions  for  general  iriflows  have  been 
obtained  by  Smith  and  Yeo  (1981)  using  methods  based  on  Hermite  polynomial  expansions. 
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In  this  paper,  we  adopt  a  different  approach  to  the  analysis  of  complex  reservoir  models. 
Based  on  a  careful  mathematical  analysis  of  die  underlying  non-linear  reservoir  model,  we  show 
that  our  watershed  model  can  be  approximated  by  a  certain  diffusion  process.  Such  diffusion 
approximations  have  been  successfully  applied  to  biology,  physics,  computer  science,  and 
many  other  areas  (see,  for  example,  Karlin  and  Taylor  (1981)).  In  particular,  such 
approximations  have  been  used  by  Yamada  (1983,  1984)  to  study  a  single  continuous  time  non¬ 
linear  finite  capacity  reservoir  with  Poisson  inflows.  In  a  similar  spirit,  Harrison  and  Shepp 
(1984)  develop  a  diffusion  approximation  for  a  cascade  of  two  discrete  time  reservoirs. 

Our  diffusion  approximation'starts  from  a  water^ed  network  of  continuous  time  nonlinear 
reservoirs  with  general  (possibly  correlated)  inflows  and  power  law  release  rules.  This  work 
differs  from  that  of  Yamada,  Harrison,  and  Shepp  both  because  we  study  an  entire  watershed 
of  reservoirs  (rather  than  a  single  reservoir)  and  because  the  power  law  release  rule  used  here 
does  not  fit  into  their  mathematical  framework.  (In  particular,  their  diffusion  limits  apply  to 
reservoirs  in  which  the  mean  inflow  rate  is  close  to  that  of  the  maximal  outflow  rate,  which  is 
assumed  to  be  finite.)  Our  approach  leads  to  a  diffusion  limit  that  is  especially  tractable  from 
a  computational  viewpoint,  primarily  because  the  approximating  diffusion  process  turns  out  to 
be  Gaussian.  In  addition,  the  mathematical  conditions  under  which  our  approximation  is  valid 
appears  to  coincide  with  a  physically  reasonable  subset  of  the  parameter  space  that  defines  the 
family  of  power  law  release  rules.  Consequently,  we  believe  that  this  approach  offers  a 
comprehensive  means  of  developing  approximations  to  brge-scale  watersheds. 

This  paper  is  organized  as  follows.  In  Section  2,  we  describe  the  basic  watershed  model, 
followed  by  a  description  of  the  diffusion  approximation  in  Section  3.  The  diffusion 
approximation  requires  that  the  user  supply  some  mean  and  covariance  information  for  the 
exogenous  inflow  processes  to  the  reservoirs  in  the  watershed.  Section  4  provides  some 
background  on  how  to  calculate  these  model  parameters.  In  Section  5,  we  discuss  the  numerical 
computation  of  the  transient  and  steady-state  distributions  of  the  approximating  diffusion  process. 
Because  the  process  is  Gaussian,  this  reduces  to  calculating  the  mean  and  covariance  of  the 
corresponding  Gaussian  random  variables.  We  show  how  these  quantities  can  be  calculated  by 
solving  a  certain  deterministic  system  of  linear  differential  equations,  the  dimension  of  which 
equals  the  number,  d,  of  reservoirs  being  studied.  Section  6  illustrates  a  different  type  of 
computation,  namely  that  associated  with  developing  an  approximation  for  the  expected  time 
required  for  a  reservoir  to  exceed  a  given  level.  Here,  one  ne^  to  solve  a  (deterministic)  linear 
partial  differential  equation  in  at  most  d  ■¥  \  variables.  In  Section  7,  we  specialize  our 
discussion  to  that  of  a  watershed  consisting  of  a  single  reservoir.  In  this  simpler  setting,  closed- 
form  formulae  are  available  for  several  quantities  of  interest  Section  8  offers  some  numerical 
evidence  which  support  the  quality  of  our  approximations,  and  Section  9  desaibes  our 
conclusions. 


2.  Description  of  the  Model 

The  drainage  network  model  will  consist  of  d  reservoirs.  In  addition,  the  following 
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notational  conventions  will  be  used: 

PQ")  -  the  collection  of  reservoirs  flowing  directly  into  reservoir  j 

R(j)  -  the  reservoirs  that  receive  inflow  water  from  reservoir  ;  ( 1  <  ;  ^  d)  (It  is  assumed 
throughout  this  paper  that  |  R(f)  |  *  IX 

Ij{t)  -  the  cumulative  exogenous  inflow  into  reservoir;  (I  sj  a  d),  and 

Sj(t)  -  the  storage  content  of  reservoir  j  at  time  t. 

We  will  assume  that  the  storage  levels  follow  "power  law"  release  rules.  The  continuity 
equation  then  states  that 


5/0  •  sj  *  tf>)  *  E  -  f „V/“J*'‘^  1  *  J  *  4  (1) 

ter(i)  J'’  J" 

where  Sj,  aj,  bj ,  Sj,  ^2*  ^3 >  positive  (deterministic)  constants.. 

We  further  assume  that  the  input  processes  are  well-behaved,  in  the  sense  that  the 
following  limits  exist  and  are  finite-valued: 


«  lim  -i-E/Zf)  (1  s  ;  s  dX  and 

!-••  t 


(2) 


C.*  -  lim  --cov(//rX  IJiO)  (1  *  j,  k  *  d)  . 

Any  input  process  with  asymptotically  stationary  increments  would  typically  satisfy  (2). 

An  example  of  such  a  network  of  reservoirs  is  given  in  Figure  1.  In  this  case  the  system 
of  equations  given  by  Equation  (1)  takes  the  following  form: 


^i(0  *  ♦  A(0  -  J[  . 

5j(r)  •  *  I^t)  -  J^*ajSj(«)‘‘du  ,  and  (3) 

5j(0  -  ^3  +  /,(0  +  J'fl,S,(tt)*'du  ♦  J['ajSj(i<)‘‘dM-  J['a,Sj(u)**du. 


3.  The  Difltasion  Approximation 

The  application  of  the  diffusion  approximation  is  quite  straightforward  and  follows  a 


sequei^oe  of  steps  which  will  be  described  in  this  section.  A  basic  assumption  of  the 
approximation  is  that  the  (1  s  i  s  (f)  appearing  in  (1)  are  small  in  a  certain  sense.  To 
calculate  the  approximations  then  requires  applying  the  following  steps: 

A-1.  Calculate  e,  where 


i/k, 

e,  =  max  a* 

Imkxd 

Our  approximation  requires  that  £ ,  should  be  "small".  Storage-runoff  models  assuming  a  small 
value  of  a  (and  hence  small  e.)  have  been  cited  in  the  literature  (see,  for  example  Laurenson, 

(1964,  p.  145)  where  the  release  rule,  r(x),  is  given  by  r(x)  =  0.0445 x*-*’  and  0.128x^-*’) 
indicating  that  this  model  may  be  a  reasonable  approximation  in  some  realistic  situations.  In  any 
case,  the  storage  continuity  equations  (1)  can  then  be  written  in  the  form 

Sy(0  -  ♦  //<)  ♦  E  C  -  [  V  w 

A  ^ 

for  1  s  ;  s  rf,  where  =  flj£.  *,  (o^  <  1,  1  s  ik  «  d). 


A-2.  Let  X.  *  e,  5^.,  Is  j%  d.  Solve  the  following  system  of  ordinary  differential  equations: 

E  Is  jsd, 

(5) 

and  such  that  x^^O)  •  Xj,  Is  js  d. 

A<3.  Solve  the  following  system  of  stochastic  differential  equations: 

^2/0  -  E  o,b,x,(^ttX(t)dt  -  ^  dBjit),  {I  sjsd)  (6) 

*€n/) 

(where  ,  B/t))  is  a  d~  dimensional  Brownian  motion  process  with  zero  drift  and 

covariance  matrix  C  «  :  I  s  i,  j  s  d)). 

We  note  that  since  (6)  comprises  a  linear  system  of  stochastic  differential  equations,  it 
follows  that  the  solution  (Zj(r),  >. ,  Zjft))  is  a  Gaussian  process.  Such  processes  possess  the 
convenient  property  that  they  are  completely  characterised  by  their  (time-dependent)  mean  and 
covariance  functions. 
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A<4.  The  approximation  is  then  given  by  the  formula 


(7) 


D 

(  •>  denotes  ^  has  approximately  the  same  distribution  as  the  process") 

It  is  clear  that  the  right  hand  side  of  (7)  inherits  the  Gaussian  structure  of  the  Zj(t)'s. 
Furthermore  the  approximating  process  appearing  in  (7)  must  necessarily  be  a  Markovian 
diffusion  process,  since  the  solution  to  the  stochastic  differential  equation  (6)  is  Markovian. 


4.  Calculating  the  Mean  and  Covariance  Structure  of  the  Exogenous  Inflow  Process 

We  now  develop  two  theorectical  models  in  which  the  mean  vector  *  (p.  :  1  ^  j  s  d) 
and  covariance  matrix  C  ^  (C.. :  1  s  i,  j  ^  d)  can  be  calculated  directly  in  terms  of  the 
underlying  "building  block"  data.  We  shall  need  to  develop  some  notation  in  order  to  precisely 
specify  our  models. 

Let  be  die  time  at  which  the  i(:*th  rainfall  event  occurs  in  the  ;’th  reservoir.  (These 
events  correspond  to  rainfall  occurring  directly  into  the  ;’th  reservoir  and  exclude  rainfall  in 
upstream  reservoirs.)  We  then  set  «  0  and  x^.  »  If  we  let  Nj(i)  be  the  number 

of  events  to  occur  by  time  t  in  the  y’th  reservoir,  it  follows  that 
Nj(t)  •  max  {n  2  0  :  7,^  s  r }.  Let  be  the  amount  of  rainfall  deposited  into  the  ;’th 
reservoir  that  propogates  into  the  reservoir;  and  we  assume  for  simplicity  that  propogation  is 
immediate.  Then  the  cumulative  exogenous  inflow  into  reservoir  j  up  to  time  t  is  given  by 


Model  1;  Assume  that  each  of  the  d  inter-event  time  sequences  :  A:  ^  1)  and  d  rainfall 
magnitude  sequences  1)  are  independent  and  identically  distributed  (i.i.d.),  and 

further  assume  that  the  2d  seqeunces  are  independent  of  one  another.  One  implication  of  this 
assumption  is  that  the  rainfall  patterns  in  each  of  the  d  subwatersheds  are  statisticdly  independent 
of  one  another.  This  is  clearly  only  reasonable  for  a  reservoir  system  covering  a  large  geographic 
area  and,  even  in  that  case,  only  in  a  very  approximate  sense.  Let 
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X.  =  l/Ex,. 

»  vzr(x,) 

v;  = 

n-  •  • 

Note  that  E[IJ(t)  I  NJ(t)]  »  Nj(tyv..  Since  Nj(t)  is  a  so-called  renewal  process,  it  follows 
that  Nj[t)  -  X-t  (See  Karlin  and  Taylor  (1975))  and  hence  as  t  -*<».  For  the 

covariances,  note  that  the  reservoir  independence  discussed  earlier  implies  that  C.^  «  0  for 
j  *  k.  As  for  Cjj,  we  use  the  variance  identity 

var(//r))  -  E[var(//0ti^<0)]  * 

(see,  for  example,  Bratley,  Fox,  and  Schrage  (1987)).  Since  |N^0)  •  »  ** 

evident  that  var(/^(/))  *  r^ENft)  +  var(Af(/)).  Again,  standard  renewal  theory  implies  that 

var(Ar<»))  -  as  r  00  (see  p.  208,  Karlin  and  Taylor  (1975)),  and  consequently. 


Cjj 


(8) 


Model  2.  In  this  model  we  will  permit  correlation  between  the  exogenous  inflows  into  each  of 
the  d  reservoirs.  In  fact,  we  will  assume  that  a  common  sequence  of  rainfall  events  affects  all 
d  reservoirs.  More  precisely,  we  will  suppose  that  for  all  1  <  /  <  d,  k  ^  1.  Of 

course,  in  this  setting,  it  is  no  longer  reasonable  to  assume  that  the  rainfall  amounts 
X^j  (1  ^  k  ^  d)  are  independent  of  one  another.  Let 


We  shall  require  that  k  *  1)  is  an  i.i.d.  sequence  independent  of 
{{X^^  :  1  %  k  s  d^  j  *  1},  which  is  itself  assumed  to  be  an  i.i.d.  sequence  (of  d-vectors).  Let 


X  ■  1/ETj 
6*  «  var(Tj) 

^ 

\  «  coyiXjp  XJ  . 

Then  the  mean  and  the  covariance  for  the  Gaussian  inflow  process  are  given  by 


*  Xv^,  and 

^jk  ’  • 
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We  conclude  this  section  with  a  brief  discussion  of  how  the  mean  vector  ^  and  the 
covariance  matrix  C  may  be  statistically  estimated  directly  from  actual  measured  exogenous 
inflow  data.  The  statistical  discussion  assumes  no  model  structure  such  as  that  associated  with 
Models  1  and  2.  (If  we  want  to  exploit  such  a  model  structure  a  more  efficient  statistical 
estimator  for  C  could  be  developed.)  Suppose  that  :  0  s  t  s  T,  I  s  j  s  d)  is  observed. 
Then,  for  1  s  y  s  d,  (l.  =  IfX)IT  is  an  appropriate  estimator  for  Py.  To  estimate  the 

covariance  matrix  C,  let  /(r)  be  the  measured  rate  at  which  the  exogenous  inflow  into  the  ;’th 
reservoir  occurs  at  time  t.  It  is  often  reasonable  to  assume  that  the  ^-dimensional  process 

((/p  ...  ,/j(0)  :  t  0)  is  stationary.  In  this  case, 

We  note  that  Cj^  is  related  to  the  cross-spectral  density  of  the  two  processes 

(/y(i)  :  r  2  0  )  and  (/^(r)  :  r  2  0  )  evaluated  at  the  frequency  0.  Q>nsequently,  estimators  from 

the  theory  of  spectral  density  estimation  can  be  used  to  calculate  for  1  s  j,  k  x  d.  (See 
Brillinger  (1981)  for  details.)  One  class  of  such  estimators  takes  the  form 

<?/.  *  IW') c6v(4(0V,{I))<i<,  (10) 

where  w^t)  is  an  appropriately  chosen  weighting  function  having  the  properties  that>Vj(r)  -*  0 
as  r  -•  oo  and  w^t)  -►  1  as  T  -•  w,  and  where 


5.  Distributional  Characterization  of  the  Approximation 

As  noted  in  Section  3,  the  diffusion  approximation  for  the  vector  storage  process 
(5(0  :  r  2  0)  is  a  Gaussian  process,  and  so  it  follows  that 

S0)°^N(m(t\K(t)) 

where  m(t)  and  K(t)  are  (respectively)  the  mean  and  the  covariance  of  the  d-dimensional  normal 
random  variable  that  approximates  S(t).  We  will  now  describe,  in  a  step-by-step  fashion,  how 
to  calculate  the  functions  m(t)  and  X(0>  thereby  completing  the  characterization  of  the 
approximation.  (For  additional  detail,  we  refer  the  reader  to  Chapter  8  of  Arnold  (1974).) 
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Let  A(t)  *  (A;y(0  :  1  ^  i,  j  s  d)ht  the  dxd  matrix  with  elements  given  by 


^(0 


iitem 
if  *  '  J.  »I><1 
0,  else, 


where  (Xj(r),  tXjit))'  is  the  solution  of  the  system  (5)  of  ordinary  differential  equations. 


The  rows  of  the  matrix  A  are  numbered  in  such  a  way  that  the  further  upstream  the 
reservoir,  the  lower  its  number;  in  other  words  we  require  that  k  6  P(f)  k  <  j  for  1  %  j  s  d. 
This  numbering  convention  will  be  used  throu^out  the  remainder  of  the  paper. 

Our  procedure  is  to: 


B<1.  Solve  the  following  deterministic  matrix-valued  differential  equation  for  <!>(/)  : 


B-2.  Calculate 


*(0  -AWfl-W 

such  that  ^0)  «  L 


rxo  -  JJ«<«)'‘C(<i<«)-‘ydu4>(ry 

where  A '  denotes  the  transpose  of  A. 

The  matrix  fXr)  can  alternatively  be  calculated  by  solving  the  matrix-valued  differential 
equation 


B-3.  Then, 


such  that  r(0)  «  0. 


in(0  -  ^£,0 

KCO  •  ijXt.t). 
c. 


(11) 


(12) 
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Of  particular  interest  is  the  steady-state  distributioD  of  the  storage.  It  is  easily  seen  that  m(r)  -•  m 
as  r  -*  00,  where  m  =  {m.  :  1  s  j  s  d)  is  the  stable  solution  of  the  differential  equation  (5)  and 
is  given  by 


(\L  *  Y, 


m.  = 


keF(j) 


a 

j 


(13) 


We  note  that  the  above  system  can  be  solved  recursively  in  closed-form  by  solving  first 
for  the  furthest  upstream  reservoirs  and  then  working  one’s  way  downstream. 

It  is  clear  that  A(t)  -*  A  as  r  -*  oo,  where 


if  te  POX 

“  j  if  A  •  j,  and 

0,  else  . 

It  is  then  reasonable  (see  (11))  to  expect  that  r(r)  -♦  P  as  r  -♦  *,  where  P  solves 

0»AP  +  PA'+C.  (1^) 

Since  A  is  lower  triangular  (due  to  the  numbering  system  adopted  at  the  beginning  of  this 
section),  and  has  negative  diagonal  elements,  it  therefore  has  strictly  negative  eigenvalues  and 
is  hence  stable.  Consequently  a  unique  solution  to  the  above  matrix  equation  is  guaranteed  (see 
Arnold  (1974)). 

If  5(<»)  is  the  (/-dimensional  random  variable  representing  the  steady-state  vector  storage 
of  the  system,  the  above  discussion  then  suggests  using  the  approximation 


S(oo)-Ar(m/E^P/Ej. 


(15) 


6.  An  Approximation  for  the  Expected  Time  for  a  Reservoir  to  Exceed  a  Given  Level 

In  this  section,  we  exploit  the  fact  that  our  approximating  process  is  a  diffusion.  (This 
is  in  contrast  to  Section  S,  where  our  analysis  relied  on  the  Gaussian  character  of  the 
approximation.)  The  Markov  structure  of  a  diffusion  process  permits  one  to  develop 
(deterministic)  partial  differential  equations  (P.D.E’s)  that  describe  a  variety  of  different 
quantitative  characteristics  of  the  process.  We  shall  illustrate  the  power  of  this  property  by 
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applying  this  idea  to  a  specific  calculation,  namely  the  development  of  an  approximation  for  the 
expected  time  for  a  reservoir  to  exceed  a  given  level.  Specifically,  we  shall  be  interested  in 
calculating  the  expected  time  required  for  reservoir  <f’s  content  to  exceed  level  p.  (We  note  that 
any  other  reservoir  can  be  viewed  as  the  furthest  downstream  reservoir  to  the  reservoirs  that  feed 
into  it.  Consequently,  if  one  is  interested  in  reservoir  j  (j  ^  d),  one  can  replace  the  original 
system  by  the  appropriate  smaller  system  of  reservoirs.  Thus,  our  assumption  is  without  loss 
of  generality.) 

If  r  is  the  time  required  for  reservoir  ds  level  to  exceed  p,  then  T  can  be  defined  as 
T  =  inf{r  2  0  :  SJt)  2  P}..  Using  the  approximation  (7),  we  conclude  that 

T  -  inf{r  2  0  :  j:j(e  ,r)/E ,  +  Z^e  j)/^  2  p} 

=  J-inf{u  2  0  :  j:»/e.  +  ZJiu)/^  2  P) 

=  JLinf{u  2  0  :  WJ^u)  2  P) 

where  W(u)  »  (Wj(u),  ...  ,Wj(u))  :  u  2  0)  is  the  diffusion  process  with  i’th  component  given  by 

Let  r,  *  inf{«  2  0  :  WJu  +  r)  2  p)  be  the  elapsed  time  for  reservoir  d’s  level  to 
exceed  p,  taken  relative  to  a  re-defined  time  origin  at  t.  Then,  we  can  set 
u{w,i)  -  E[r,  I  W(r)  *  w]  for  r  2  0,  w  6  E^.  The  diffusion  W(')  has  infinitesimal  generator 
L  given  by  the  linear  partial  differential  operator 

L  *  ^  1  ^  ^ij 

..I  £ .  *0-0)  ^  ^ 

The  theory  of  diffusion  processes  (see  Karlin  and  Taylor  (1981))  establishes  that  u 
satisfies  the  P.D.E.  given  by 

♦  L«  =  -1  (16) 

dt 
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subject  to  u(w,t)  =  0  for  a  p.  Hence,  we  can  calculate  an  approximation  to  the  expected 
time  E(2)  required  for  reservoir  d’s  level  to  exceed  p  by  following  the  procedure: 

C-1.  Solve  the  (deterministic)  P.D.E 


du 

It 


f: 


du 

f:  \ 

dw. 

+  J.  V 

2^  E,  dw.dw. 


=  -1 


for  u(w,t)  in  the  domain  w  6  H'',  r  s  0,  subject  to  M(>v,f)  =  0  for  fe  p. 


C-2.  Then  E(7)  -  0). 

E . 


(17) 


7.  The  Single  Reservoir  Case. 

The  mass  balance  equation  for  a  single  reservoir  is  given  by  the  following  relation 

S(0  -  +  7(0  -  J['flS(u)‘du 

where  Sq,  a,  b  >  0,  and  7(r)  is  a  non-decreasing  process  that  represents  the  cumulative  input 
to  time  t.  As  required  by  assumption  (2),  we  assume  that  the  limits  ^  -  lim  EI{t)/t  and 

C  »  lim  var/(r)/r  exist. 

Following  the  procedure  outlined  in  Section  3,  it  is  first  necessary  to  solve  the  ordinary 
differential  equation  given  by 


x(r)  *  n  -  (xx(r)*,  with 

where  =  a  ^  and  a  =  1.  The  solution  to  the  di&erential  equation  is  then  given  by 

where 


9(')  =  I 


du 
^  -  u 


7  ■ 
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In  particular,  if  6  =1,  then 


/. 


(p(x)  -  lod 


1*  “-*0 


\ 


and  so  the  solution  of  the  differential  equation  becomes  x(t)  •  \i  x^.  The  next  step 

in  the  procedure  is  to  solve  the  stochastic  differential  equation  for  the  diffusion  approximation 

which  is  given  by  dZ(t)  *  -6x(r)*‘^Z(/)<ir  +  dB(t)  where  S(‘)  is  a  Brownian  motion  having  a 
variance  parameter  C.  Then,  in  the  linear  case,  we  can  apply  (11)  to  obtain  the  result 
r(r)  «  C(1  -  exp(-2r))/2.  Since  e.  *  a,  the  transient  approximation  (7)  becomes 


S(0  -  N((vi  -  e-^di  -  s^))/a,  C(1  -  e'^) 


(18) 


In  the  case  where  b  is  not  equal  to  one,  we  can  easily  supply  a  closed-fonn  solution  for 

the  steady-state  storage  S(«).  In  particular,  (14)  takes  the  form  T  ■  Then  (15) 

yields 


(19) 


8.  Numerical  Examines 

We  now  consider  two  numerical  examples.  In  the  first  case  we  will  compare  the  diffusion 
approximation  to  exact  known  results  from  Harristm  and  Resnick  (1976),  and  in  the  second  case 
we  will  compare  the  diffusion  approximatitm  to  simulations  for  a  cascade  of  two  reservoirs. 

Example  1.  The  behaviour  of  the  system  for  a  single  reservoir  with  Poisson  inflows  and 
exponentially  distributed  jump  sizes  will  first  be  examined.  The  statistics  of  the  inflow  process 
can  then  be  calculated  as  in  Section  (4).  This  case  has  been  solved  in  closed  form  by  Harrison 
and  Resnick  in  the  ^>ecial  case  where  the  storage  system  empties  in  finite  time  i.e.  for  the  case 
where  b  <  1.  In  this  case  define  two  functions  p(x)  and  K{x)  by  the  formulae 

p(x)  ■  X/(var  *),  and 

K(x)  -  vp(x)exp(-vj['(l  -p(/))dt)  , 


where  X.  and  v  are  defined  as  in  Model  1  of  Section  4. 
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The  steady 'State  density,  fj(x),  of  the  storage  process  is  then  given  by  the  formula 


fjix)  .  /C(ry(l  *  ^’K(m  .  (20) 

Similarly  the  density  of  the  runoff  process  fj^x)  is  given  by  the  iV  mula 

fj[x)  «  abx^~^fj[ax^).  It  is  clear  from  th^  formulae  and  from  the  fact  that  the  storage  is 
depleted  in  finite  time  that  the  distributions  of  the  storage  and  runoff  processes  have  atoms  at  0. 
A  comparison  of  the  densities  of  storage  for  the  diffusion  process  and  for  the  exact  solution  of 
the  storage  process  for  various  values  of  a  is  shown  in  Figure  2. 

Example  2.  Consider  a  cascade  of  two  storage  systems  fed  by  a  single  exogenous  inflow  process 
/^(*)into  the  first  storage  system.  In  this  case  the  matrix  A  (see  Section  5)  takes  the  form 


A 


k.-l 

1 


0 


(21) 


We  can  calculate  the  steady-state  variance  of  the  storage  process  by  solving  Equation  (14) 
for  r,  which  results  in  the  following  formula 


'  1 

1 

‘Cn 

^1 

2 

1 

^11 

Ajj**"  Ajj 

(22) 


where  «  lim  w  I  and  the  A..'s  are  the  corresponding  elements  of  the  matrix  A  defined 

in  (21).  The  steady-state  distribution  of  this  non-linear  cascade  can  then  be  approximated  as  in 
(15). 


The  diffusion  approximation  for  the  cascade  was  compared  with  simulation  runs  for  the 
process.  The  simulation  model  was  based  on  approximating  the  continuous  time  model  given 
by  equation  (1)  with  the  following  approximate  discrete  time  model 

Si(n+1)  -  5i(n)  -  /,(Af)  -  ^23) 

Sj(n+1)  -  S^n)  ■  ajAtS^(n*lf'  -  . 

where  {5j(/:)  :  n  2  0}  and  {Sjin)  :  /i  2  0}  are  the  storages  of  the  first  and  second  system 
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J 


respectively.  Such  a  discretization  is  necessary  in  order  to  solve  the  system  of  integral  equations 
(1)  numerically.  The  inflow  process  :  n  2  0}  was  simulated  using  the  formula 

AT/A.) 

/.(A,)  .  5;  C.,  (M) 

>1 

where  the  G^  .  are  independent  and  identically  distributed  Gamma  random  variables  and  the 
are  independent  and  identically  distributed  Poisson  random  variables  with  mean  XAr. 

In  our  cascade  example,  the  parameter  X  for  the  Poisson  process  :  j  2  0}  was 

assumed  to  be  one,  as  was  the  mean,  v ,  and  the  standard  deviation,  t],  of  the  gamma  deviates, 
N^j.  In  this  case,  is  equal  to  two,  and  Cj^,  and  Cjj  afc  all  zero.  The  time 

parameter  A/  was  put  equal  to  0.1  and  the  storage  process  simulation  was  run  for  10000  time 
steps  (to  t  s  1000).  Examples  of  the  resulting  sample  paths  for  runoff  from  the  two  storage 
.systems  are  shown  in  Figures  3  and  4,  and  the  results  of  simulations  are  provided  in  Hgure  5. 
Confidence  intervals  for  the  estimated  parameters  were  calculated  using  an  assumed  t  distribution 
based  on  ten  replications  of  each  simulation.  The  justification  for  this  confidence  interval 
methodology  is  given  in  Law  and  Kelton  (1981). 


9.  Conclusions 

In  this  paper,  we  have  developed  a  diffusion  approximation  for  a  watershed  consisting 
of  d  reservoirs  with  nonlinear  power  law  release  rules.  Our  simulation  experiments  suggest  that 
it  works  reasonably  well  for  some  values  of  the  release  rule  parameters  which  are  physically 
plausible.  Of  course,  the  assumptions  guaranteeing  the  validity  of  our  approximation  are  not 
universally  applicable,  and  further  research  is  therefore  necessary  to  fully  delineate  its  domain 
of  applicability. 
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Appendix  •  Justification  for  the  Diffusion  Approximation 

We  wish  to  show  that  when  the  parameter  e,  is  small,  then  the  diffusion  approximation 
is  reasonable.  To  formulate  this  problem  mathematically,  let 

5/e,  0  *  5/e,  0)  +  7/0  +  53  f 'c^£**S/E,  uf^du  -  f'cLe^'^/E,  uf'du,  1  s  ;  s  d, 

*eflw  J'’ 

(25) 

be  the  storage  contents  processes  associated  with  the  parameter  e .  We  need  to  show  that 
S/e,  t)  can  be  approximated  as  in  (A>4)  when  e  i  0. 

To  study  the  behaviour  of  {S/e,  t)  :  1  s  j  ^  t  ^  0}  as  t  |  0,  we  need  to  strengthen 
the  hypothesis  (2)  on  the  inflow  processes  somewhat.  (In  particular,  we  need  to  say  something 
about  the  distribution  of  the  inflow  processes.)  In  addition,  since  our  approximation  requires 
a  re-scaling  of  the  spatial  co-ordinates,  our  initial  condition  for  S/e,  0)  ought  to  reflect  this. 
Our  precise  mathematical  requirement  is: 

e*‘^(eS/e,0)  -  Xj,  ,  eS/e,0)  -  x^,  Ef/r/E)  -  Pj/,  -  ,  eZ/z/e)  -  p/] 

(Z/0), -.  ,  Z/0),  B/0.  -  .  5/0], 

(26) 

as  E  I  0,  where  denotes  convergence  in  distribution  and  the  process  (B/rX***,  B/r))  is 
assumed  to  be  a  Brownian  motion  with  zero  mean  and  covariance  matrix  C.  From  a  practical 
standpoint,  the  stengthening  of  (2)  to  the  above  distributional  requirement  is  quite  mild. 
Virtually  all  stochastic  processes  that  exhibit  the  behaviour  (2)  also  exhibit  the  Brownian 
behaviour  assumed  above  (Ethier  and  Kurtz  (1986)). 

Under  the  assumption  (26),  we  can  now  proceed  to  analyze  S/e,  *)asE  |  0.  Note  that 
(25)  can  be  re-written  as 

ES/E,r/E)  »  eS/e,0)  ♦  e//r/E)  ♦  5}  |V^5/e,«/E))**du  -  ^'aftSft,u/t)f*du  . 

By  assumption,  we  have  eZ/r/e)  p^r  as  e  i  0,  and  so  assuming  in  addition  that  the  limit 
process 


x/r)  «  lime  S/e, r/E) 

«i0  ' 


(27) 
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exists  for  I  x  j  x  d,  the  limit  must  satisfy  the  integral  equation 

^;<0  *  X.  +  M/  +  E  [  1  s  ;  s  rf. 

Expressed  as  a  differential  equation,  the  abow  result  becomes 

“  1*;  +  E  ViCO*"  -  and 

teJTi) 

Xj(0)‘Xj,  l^J^d. 

In  order  to  establish  the  existence  of  the  limit,  one  needs  to  prove  the  relative 
compactness  of  the  ^ily  of  processes  {e  Sj(t,  t/t)  :  e  k  0).  This  can  be  done  by  appealing 
to  standard  tools  from,  for  example,  Ethier  and  Kurtz  (1986).  We  can  view  (27)  as  a  law  of 
large  numbers  for  Sft,  *)•  develop  a  "Central  Limit  Theorem"  we  re-write  (25)  as 

e-“(tS/t,i/0 -xffi]  .  e-“[eS,<t.O) -x)  - 

*  E  f -  xfufyu 

law* 

-  -  x0)iu. 

A  Taylor  series  expansion  to  one  term  then  establishes  that 

e-‘«((eS^t,Wt))*‘  -x/uy-l  -x,(«))] 

where  ^|(e,u)  lies  between  xj^u)  and  e  and  hence  tends  to  Xj^(u)  as  e  i  0. 

Assume,  for  the  moment,  that  we  can  establish  the  existence  of  the  limit 

-  Xj(r),~,eSj(c,t/e)  -  xj(0)  (Zj(a~^j(0)  • 

Since  e'‘'*[eS,(e,0)  -  Xj,~,eSj(e,0)  -  xj [Zj(0),~,Zj(0)l,  the  limit  process 
(Zj(/X  -  must  then  satisfy  the  following  stochastic  differential  equation 

z/o  ■  Z/O)  *  Bfi)  .  £  ^'aJ,fJiut"Ziu)<lu  -  ^'afiicf.uf'-'Zfuyiu  . 

As  before,  the  existence  of  the  limit  may  be  justified  by  using  a  relative  cornpacmess 
argument  (Ethier  and  Kurtz  (1986)). 
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Figure  1.  A  storage  system  consisting  of 
three  reservoirs  where  reservoirs  one  and 
two  feed  into  reservoir  three. 


Figure  2.  A  Comparison  of  the  exact 
solution  and  the  diffusion  approximation  for 
a  single  reservoir  (b  =  0.5,  X  =  1.0,  v  =  1.0) 


Figure  3.  Sample  path  of  first  reservoir  fed  by 
compound  Poisson  inflows  (rate  1.0)  with 
gamma  distributed  jumps  ( mean  jump  size  2.0) 
and  with  Ar  ^  0.1,  U]  s  0.5,  b,  «  1.5. 


Figure  4.  Sample  path  of  second  reservoir  with 
same  parameters  as  first  reservoir  (see  Figure 


Figure  5.  Exceedence  probabilities  for 
runoff  from  second  reservoir  with  Uj  s  s 
0.3;  Solid  line  represents  diffusion  and  error 
bars  represent  95%  CJ.  for  simulation. 
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